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Abstract 

An aerodynamic design algorithm for turbulent nows using unstructured grids is described. The current approach uses 
adjoint (costate) variables to obtain derivatives of the cost function. The solution of the adjoint equations is obtained by 
using an implicit formulation in which the turbulence model is fully coupled with the flow equations when solving for the 
costate variables. The accuracy of the derivatives is demonstrated by comparison with finite-difference gradients and a few 
sample computations are shown. In addition, a user interface is described that significantly reduces the time required to set 
up the design problems. Recommendations on directions of further research into the Navier-Stokes design process are 
made. 


Introduction 

Because of rapid advances in computer speeds, and improvements in 
flow-solver and grid-generation algorithms, a renewed emphasis has 
been placed on extending computational fluid dynamics (CFD) beyond 
its traditional role as an analysis tool to design optimization. Among the 
methodologies often employed are gradient-based techniques, in which a 
specified objective is minimized. In this framework, the gradients of the 
objective function with respect to the design variables are used to update 
the design variables in a systematic manner to reduce the cost function 
and to arrive at a local minimum. Many techniques have been used to 
obtain the necessary derivatives, including finite differences, direct dif- 
ferentiation, and adjoint methods. Many of the methodologies and im- 
plementations are discussed in Refs. 1, 4, 7, 8, 11-13, 17, 19-24, 29, 30 
and 34. 

Although most of the above mentioned references deal with inviscid 
flows, a few have addressed viscous computations of turbulent flows. In 
Ref. 18, Hou et al. used a direct differentiation approach in which the de- 
rivatives of the flow solver were obtained with ADIFOR. 9 In Ref. 18, the 
turbulence model used was the Baldwin-Lomax 6 algebraic model, which 
was differentiated along with the flow equations. Jameson recently de- 
veloped a design methodology for turbulent flows based on an adjoint 
formulation. 21 Here, the Baldwin-Lomax turbulence model was also em- 
ployed but was assumed constant and was therefore not differentiated. 
This same assumption was also recently used in the work of Soemar- 
woto. 34 

For unstructured grids, the work to date has been primarily focused on 
inviscid computations in both two and three dimensions. 1,8,13,14 ’ 27,28 In 
Ref. 1, the adjoint equations and boundary conditions were derived for 
the incompressible Navier-Stokes equations, and some design examples 
were demonstrated. However, turbulence effects were not included. In 
the work of Mohammadi, 26 two-dimensional Navier-Stokes results were 
presented in which turbulence effects were included. In this reference, 
automatic differentiation was used to differentiate the necessary compo- 
nents of the flow solver. 

The purpose of the present study is to extend the work in Ref. 1 to the 
compressible Navier-Stokes equations, including a fully coupled field- 
equation turbulence model. However, because a continuous adjoint ap- 
proach for unstructured grids requires accurately computed second de- 
rivatives, 1 a discrete adjoint approach is used in the present study. The 
methodology is discussed, and the accuracy of the derivatives is estab- 
lished. A few design examples are given to demonstrate the technology. 
Also, a user interface has been developed to facilitate setup of the de- 
signs, and a description of the interface is included. 


Nomenclature 


A 

a 

C* 

C] 


area of control volume 
speed of sound 

constant used in Sutherland’s law for viscosity 
lift coefficient 


Cj drag coefficient 

c , Ck c constants used in Spalart-Allmaras turbulence 

V b z’ v i 

model 

c c c constants used in Spalart-Allmaras turbulence 

w ] ’ w 3 

mode 


D vector of design variables 

Dj component i of design variable vector 


d 

E 


distance to nearest surface 

total energy per unit volume 

fluxes of mass, momentum, and energy 


inviscid contribution to fluxes 
f v viscous contribution to fluxes 


f, g components of inviscid fluxes 

f v , g v components of viscous fluxes 



I 

Ic 


K 


functions used in the turbulence model 

functions used in the turbulence model 

augmented cost to be minimized 
cost to be minimized 

Karman constant 


M„ 

N 

n 

Pr 

Pr, 


free-stream Mach number 

B -spline basis functions 

unit normal to boundary of control volume 

Prandtl number 

turbulent Prandtl number 
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p pressure 

Q vector of dependent variables 

q x , q y components of heat flux 

R residual for a control volume 

Re Reynolds number 

S magnitude of vorticity 

s parameterization variable for B-splines 

T temperature 

t time 


A 5t* + A<H = 0 

dn do 


y + wall coordinate 

a angle of attack 

dO boundary of control volume 

y ratio of specific heats 

^ laminar viscosity 

|X t turbulent viscosity 

x> |i/p 

u, ii,/p 


(i) 


where A is the outward-pointing normal to the control volume. The 
vector of dependent variables Q and the flux vectors ft and ft are 
given as 


Q = 


( 2 ) 


u 

magnitude of velocity 


pu 


pv 

u, v 

X 

Cartesian components of velocity 
grid-point locations 

ft = fi+gj = 

pu 2 + p 
puv 

i + 

pvu 
pv 2 + p 

x,y 

Cartesian coordinates 


(E + p)u 


(E + p)v 


(3) 


and 


fv«+g,j = 


0 


0 

x 


x 


i + 

*y 

X 


x 



yy 

UT„ + vx xy - q x 


uT yx + vT yy -q y 


j (4) 


Here, ft and F v are the inviscid and viscous flux vectors respec- 
tively; the shear stress and heat conduction terms are given as 


o 

dependent variable for turbulence model 

M„2 

^ = (H + H,)j^3(2u x -v y ) 

(5) 

p 

density 



G 

XXX 
‘xx’ fc xy» v yy 

constant for turbulence model 
shear stress terms 

M„2 

T yy = (p + p,)— -(2v y -u x ) 

(6) 

¥ 

costate variables 

= V = (l l + l i .)R^( u y + v «) 

(7) 

Superscripts: 

dimensional quantity 

fp p, W 

q * Re(y- l)lPr + Prjax 

(8) 

Subscripts: 

variation 

^ p, ^a 2 

Qy " Re(y - 1 )lPr + Prjdy 

(9) 

DO 

free-stream quantities 

The equations are closed with the equation of state for a perfect gas 


Governing Equations 

, i\fr: (U 2 + V 2 )l 

p = (Y 1)1 E p' 2 

(10) 


Flow equations 

The governing equations are the time-dependent Reynolds-averaged 
Navier-Stokes equations. The equations are expressed as a system of 
conservation laws that relate the time rate of change of mass, momen- 
tum, and energy in a control volume of area A to the spatial fluxes of 
these quantities through the volume. The equations (nondimensional- 
ized by free-stream density, speed of sound, temperature, viscosity, 
thermal conductivity, and a reference length) are given as 


and the laminar viscosity is determined through Sutherland’s law: 


jl _ (1 +C*)(f/f^) 3/; 

A- - + 


(ii) 


where C* = 198.6/460.0 is Sutherland’s constant divided by a 
free-stream reference temperature, which is assumed to be 460° R. 
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Turbulence model 

For the current study, the turbulence model of Spalart-Allmaras is 
used. 35 This is a one -equation turbulence model given as 


^ = ^{V.[( V + (l + c b 2 )X))V{,]-c b iV^} 

- - 5 f '.)(s/ * S< 1 - V 5 " + S i . , 'i au ’ 


(12) 


ables). In Eq. (21), I C (Q, D) represents the cost that is to be mini- 
mized. Examples of suitable cost functions include the difference 
between the lift coefficient for the airfoil and a desired lift, the drag 
coefficient, and the difference between the pressure distribution and 
a desired pressure distribution. 

The variation of Eq. (21) is given by 


61 = 


3I C - 3I C - T 

< Q + _LD + 'F 
3Q V 3D 


r3R - 

l5Q Q * 


(9r aRax^i 

Ud + 3X3DJ . 


( 22 ) 


where 


f 


X 3 

X 3 + c v , 


The terms involving Q can be eliminated by regrouping terms 
and requiring the coefficients of Q to vanish; the costate variables 
(13) are the solution of the following equations 


X = 


\) 

v 


(14) 



(23) 


S 


S + 


M, j> 

Re K 2 d 2 v 2 


The remaining terms for the variation in the cost function are then 
(15) given by 


and 


61 


31 C , ,,, T f3R 3R9XY) 

,3d Ud 3X3dJJ 


D 


(24) 


f 


= i 2_ 


(16) 


In these equations, S is the magnitude of the vorticity, and d is the 
distance to the nearest wall. The function f w is given as 


= g 


1 + c£. 'l 


■V 


(17) 


where 


g = r + c Wj (r 6 -r) 


and 


M r _6_ 

Rested 2 


(18) 


(19) 


The last term in Eq. (12) is used when specifying the transition loca- 
tion. Although the flow solver includes this term, the computations in 
the present paper are all assumed to be fully turbulent, so this term is 
not used. Therefore, the definitions of f t and f 2 , which are associated 
with these terms, are not given. After kq. (12) is solved for “6, the 
eddy viscosity is computed as 

n, = pu, = p\)f V| (20) 

Adjoint Equations 

In the adjoint approach for design optimization, a cost function is 
defined and augmented with the flow equations as constraints: 

l[Q,D,'F,X(D)] = 1 C (Q, D) + l P T R[Q, D, X(D)] (21) 

where R represents the vector of discrete residuals, X is the location 
of the grid points, D is the vector of design variables, and T are the 
Lagrange multipliers (also referred to as the costate or adjoint vari- 


After the costate variables are determined from Eq. (23), they are 
used in Eq. (24) to obtain the sensitivity derivatives. Note that this 
process requires the solution of both the flow equations and the cos- 
tate equations. However, the derivatives of the cost function with re- 
spect to all design variables are obtained independently of the num- 
ber of design variables. 

By examining Eqs. (5)-(9) along with Eqs. (12)-(20), it is apparent 
that the solution of the flow equations and the turbulent viscosity are 
highly dependent on one another. Therefore, the vector of residuals 
that require linearization in Eqs. (23) and (24) includes the contribu- 
tions from both the flow equations and the turbulence model. Like- 
wise, the dependent variables, Q , include the conserved flow vari- 
ables as well as v so that solving for the costate variables with Eq. 
(23) requires the solution of a block 5x5 system of equations for 
two-dimensional calculations and a 6x6 system in three-dimensions. 

Many of the terms in Eqs. (12)-(20) have a complex dependency 
on the dependent variables, the design variables, and the distance to 
the wall; these terms must be accurately differentiated in order to ob- 
tain accurate derivatives. In the present work, the differentiation of 
both the flow equations and the turbulence model is accomplished 
by “hand differentiating” the code. Although this procedure is some- 
what tedious, experiments in which the eddy viscosity was assumed 
to be constant (and, therefore, not differentiated) yielded very poor 
accuracy with many derivatives of incorrect sign when compared 
with gradients obtained with finite differences. The strong coupling 
of the flow equations and the turbulence model is in contrast to Refs. 
21 and 34 where the constant viscosity assumption was used. How- 
ever, in those references, an algebraic turbulence model is used, 
whereas here, a field equation is solved to obtain the eddy viscosity. 

In Eq. (24), the terms that involve *F T [( 3R/3X ) ■ (3X/3D)]D 
represent the change in the cost function that results from a change 
in the grid. Reference 1 shows that the contributions of these terms 
diminish as the grid is refined except at geometric singularities such 
as trailing edges. Because the position of the trailing edge is fixed in 
the present work, these terms are currently not included in the com- 
putations so that Eq. (24) can be evaluated by looping over only a 
small subset of edges in the mesh rather than the entire mesh. How- 
ever, in order to effectively deal with cases that require derivatives at 
the trailing edge and to better ensure convergence of the optimiza- 
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tion procedure, even on coarse grids, it is recommended that these 
terms be included. 


Solution Procedures 

For the flow equations, the inviscid flux contributions are evaluated 
by using an approximate Riemann solver, 31 and the viscous contribu- 
tions are discretized with a central -difference approach. The solution is 
obtained by using an implicit solution methodology with multigrid ac- 
celeration. Details may be found in Refs. 2, 3, and 10. The adjoint 
equations are a linear system of equations that can be solved with a 
technique such as preconditioned GMRES. 32 However, in this work, a 
time derivative is added to the equations so that they can be solved 
with a time-marching procedure. The motivation for adding the time 
term is that this approach often converges in situations for which the 
preconditioned GMRES might otherwise fail. This feature is particu- 
larly useful when the turbulence model is fully coupled because the 
turbulence production term tends to reduce diagonal dominance. Be- 
cause the adjoint equations represent a linear system of equations, the 
matrix-vector products are currently formed by simply passing the 
vector to the residual routine in place of the costate variables. By form- 
ing the matrix-vector products in this way, the largest contribution to 
memory requirements is through the preconditioner (incomplete 
lower/upper (LU) decomposition with no fill, ILU(O)), so that the re- 
sulting scheme requires roughly the same amount of memory as the 
flow solver. Note that this procedure essentially requires recomputa- 
tion of the linearization of the residual for each matrix-vector product. 


is not drastically effected by the choice of step size. For this grid 
under these flow conditions, the maximum turbulent viscosity in the 
flow field is approximately 2500. 



Figure 1 . Grid used for studying accuracy of derivatives. 


Grid Generation and Mesh Movement 

The unstructured meshes used in this work are generated with the 
software package described in Ref. 25. This employs an advancing 
front method that generates good quality grids for both inviscid and 
viscous calculations. 

During the design process, the mesh is continuously updated as the 
shape of the geometry changes. This is accomplished using the tech- 
nique described in Ref. 1 , which shifts nodes near viscous surfaces by 
interpolating the changes in the coordinates at the end points of the 
nearest surface edge. This technique is blended with a smoothing pro- 
cedure so that away from the highly stretched cells near the surface the 
mesh movement reverts to that of the smoothing/edge- swapping pro- 
cedure described in Ref. 38. The combined procedure has been found 
to work well for viscous grids with highly stretched triangles and very 
close spacing normal to the wall. Further details can be found in Ref. 
1 . 


Accuracy of Derivatives 

The accuracy of the derivatives is established by comparing results 
obtained by using the adjoint formulation with finite-difference deriv- 
atives. The case considered here is a symmetric airfoil at a free-stream 
Mach number of 0.55, an angle of attack of 1° , and a Reynolds num- 
ber of 9 million, based on the chord of the airfoil. The grid used is 
fairly coarse with only 2700 nodes and a spacing at the wall of about 
5 x 10^ . (See Fig. 1 .) The spacing at the wall has been chosen to be 
large enough so that a stretched mesh can be obtained while allowing 
the surface to be perturbed without moving the interior grid points. 
This is done in order to remain consistent with the assumption that the 
interior mesh sensitivities are neglected. Although this grid is obvi- 
ously inadequate for resolving the boundary layer accurately, it is suf- 
ficient for verifying the linearization of the flow solver. When the gra- 
dients are computed with finite differences, a central -difference 
formula is used with a fixed step size for each design variable, and all 
computations are converged to machine accuracy. For grids in which 
closer spacing at the wall is used, Hou et al. 18 have shown that obtain- 
ing derivatives from finite differences can be highly sensitive to the 
step size and to the level of convergence of the flow solver. With the 
spacing at the wall used here, the flow solver is easily converged to 
machine zero, and numerical experiments indicate that the derivative 


For this test, the geometry of the airfoil is described with a third 
order B -spline with 39 control points. The derivatives of lift and 
drag coefficients with respect to angle of attack and to three shape 
design variables are evaluated. The shape design variables corre- 
spond to the y-position of three control points located at 
x/c = 0.103, x/c = 0.789 , and x/c = 0.972 and designated as 
Dj , D 2 , and D 3 , respectively. As seen in the tables below, the de- 
rivatives obtained with the adjoint approach are in good agreement 
with the finite-difference derivatives. 


Table 1. Accuracy of Derivatives for Lift Coefficient 



Finite-difference 

Adjoint 

dc l 

da 

5.8278 

5.8278 

3c, 

3d; 

-1.9065 

-1.9067 

3c, 

3D^ 

1.3505 

1.3505 

3c, 

35 ; 

0.44746 

0.44746 
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Table 2. Accuracy of Derivatives for Drag Coefficient 



Finite-difference 

Adjoint 

9c d 

9a 

0.057503 

0.057503 

9c d 

9D, 

-0.50394 

-0.50401 

9c d 

9D 2 

-0.063550 

-0.063547 

9c d 

9D, 

-0.0084115 

-0.0084114 


Surface Representation and Graphical Interface 

In the current study, the geometries are modeled with B-splines, 
which offer great flexibility in the definition of the surfaces. By vary- 
ing the polynomial degree and the number of control points, a wide 
range in the number of design variables and in surface fidelity can be 
obtained. On one hand, the design variables can be made to correspond 
to the individual grid points on the surface by choosing a linear poly- 
nomial and an appropriate number of control points. Conversely, a sin- 
gle polynomial curve of degree M (known as a Bezier curve) can be 
used to describe the geometry by choosing the number of control 
points to be M + 1 . In addition, through the knot sequence associated 
with the spline, curves with sharp breaks in the surface, such as those 
that occur in cove regions and blunt trailing edges, can still be repre- 
sented in a single curve. 

Spline fitting of input coordinates 

Rather than using a conventional cubic spline of the input coordi- 
nates, a B-spline of specified order and with a specified number of 
control points is matched to the input coordinates with a least-squares 
procedure. The design variables are, then, the coordinates of the B- 
spline control points, which can be considerably fewer in number and 
are more geometrically meaningful than the original input coordinates. 
The following is a description of the B-spline representation and the 
least-squares procedure. 

B-spline curves are described in detail in Ref. 15. They are defined 
as the sum of products of con U*ol -point coordinates and corresponding 
basis functions. The basis functions depend on the parametrization of 
the spline and a knot sequence and are defined recursively as follows: 


N j ( S ) : 


S - S: 


. N l\ s) + huL^ N: ~i (s) 


j 1, s,_, <s<Si 
1 0, else 


(25) 


where n is the degree of the basis function. The minimum and the 
maximum values of the parameter s appear n times at the beginning 
and the end of the knot sequence, respectively, so that the first and last 
control points correspond to the end points of the B-spline. 

A uniform parameterization is formed by setting the parameter s 
that corresponds to each input coordinate equal to the number of the 
coordinate in the sequence, starting from zero: 

sj = j, je [0, M'l (26) 

where M' is the number of input coordinates. The knot sequence is 
formed by uniform division of the parameter space. At each of the Sj , 


each of the M basis functions N" is computed, which forms an 
M'xM matrix. 

Given the values of the basis functions at each input coordinate, 
an overdetermined linear system is obtained: 

M 

Xj = ^XiN^Sj) (27) 

i = 0 

where xj is the y-th input coordinate and M is the number of control 
points. The \ } are the unknown control point coordinates, and the 
xj are the input coordinates. The first and last control points are set 
equal to the first and last input coordinates, and the corresponding 
equations are removed from the system. The system is then solved in 
the least-squares sense by using Householder transformations as de- 
scribed in Ref. 16. After fitting each segment of a curve, the B-spline 
segments are concatenated into a single B-spline by concatenating 
the knot sequences and merging the control point coordinates. 

Graphical interface 

A user interface has been written to facilitate the setup of each de- 
sign. The intent is to not only speed up the process of setting up each 
case but to help eliminate errors. Because a B-spline representation 
is used to describe the geometries, the points on the surface of the 
airfoil must lie on this surface. Therefore, a least-squares fit of a B- 
spline curve to the point description of the airfoil is used. This pro- 
cess is described in a previous section and is demonstrated in Fig. 2 
for a two-element airfoil, where the positions of the control points 
are indicated by the symbols. The control panel on the right shows 
that for the second element (the flap) a third-order B-spline with 25 
control points is used to define the surface. The positions of these 
control points are subsequently used as the design variables. Note 
that although a cubic representation is used in the present example 
the order of the spline can range from linear to M - 1 , where M is 
the number of points that define the airfoil. In the case for which M 
points are used to define a curve of degree M - 1 , the resulting 
curve corresponds to a Bezier representation. In the present exam 
pie, the B-spline is evaluated at 129 points, which are then used as 
input to the grid generation to define the surface. Although not 
shown, a similar procedure has also been used for the main element. 
In Fig. 2, a “tear off' menu is also shown, which allows the user to 
choose whether various symbols are displayed. 

In the next figure (Fig. 3), upper and lower limits have been 
placed on the y-coordinate of many of the control points, with the 
limits depicted by the extent of the lines above and below the current 
placement of the design variable. The limits are used as side con- 
straints during the optimization process to prevent the occurrence of 
nonphysical geometries during the design process. The x-coordi- 
nates can also be chosen as design variables by placing distinct 
upper and lower bounds on their positions; if the upper and lower 
bounds are the same, then the design variable is not allowed to 
change and is not considered in the optimization process. 

Other capabilities include the ability to zoom and to translate the 
geometry to manipulate it into position. This capability is particu- 
larly useful when placing side constraints near the trailing edges of 
the airfoils. Also, the shape of the airfoil surfaces can be changed 
simply by picking and moving the control point, and the initial posi- 
tion can be recovered. 

After splining the surface and setting the limits of the design vari- 
ables, output is written in a format that is suitable for the grid gener- 
ation program. A file that contains the geometric information, such 
as the positions and limits of the control points, is also written. This 
file is continually updated during the design process to reflect the 
changing shape. This file can be read in at a later date to reset the 
limits on the design variables, to add more points to the surface defi- 
nition, or to reshape the geometry. 
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Figure 3. Setting the upper and lower limits on design variables. 
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Optimizer 

The optimizer used in the current study is KSOPT, 39 which uses a 
quasi-Newton method to determine the search directions and a polyno- 
mial line search technique to determine the step length in the descent 
direction. This code was chosen because it is capable of multipoint de- 
sign and can handle both equality and inequality constraints. In addi- 
tion, upper and lower bounds can be placed on design variables; this 
approach is currently used to enforce the geometric constraints neces- 
sary to maintain a viable geometry throughout the design cycle. 

Results 

Two sample results are given below. The first case is a computation 
of the flow over an airfoil at a free -stream Mach number of 0.4, an 
angle of attack of 2° , and a Reynolds number of 5 million. The goal of 
the computation is simply to obtain a specified pressure distribution. 
The grid used for this computation is shown in Fig. 4 and consists of 
approximately 5500 nodes with 128 nodes on the surface of the airfoil. 
The spacing at the wall is 1 x 10~ 5 of the chord length yielding a y + 
of about 2. For this case, a single eighth-order Bezier curve is used to 
parameterize the surface, and only three design variables are allowed 
to change during the design process. The geometry is perturbed by dis- 
placing three of the control points in the initial B-spline definition, and 
the solution over this geometry is used for the target pressures. After 
10 design cycles, the cost function is reduced from approximately 
1.5 x 10~ l to 3.0 x 10 7 , and the root mean square of the gradients is 
reduced from 1.4 to 5.8 x 10^ . The initial and final pressure distribu- 
tions and geometries are shown in Figs. 5 and 6. As seen, the target 
pressure distribution is obtained, and the geometry returns to that of 
the airfoil in the perturbed position. 



Figure 4. Grid used for first test case. 



x/c 

Figure 5. Initial and final pressure distributions for case 1 . 



x/c 

Figure 6. Initial and final geometries for case 1 . 


The next case is also an inverse design case in which the objective 
is to match a specified pressure distribution. However, this case is 
significantly more difficult because 75 design variables are used. 
The initial grid used for this case consists of about 5300 nodes, and 
has a spacing normal to the wall of 1.0 x 10~ 5 . (See Fig. 7.) The 
free-stream Mach number is 0.725, the angle of attack is 2.54° , and 
the Reynolds number is 6.5 million. For this case, 20 design cycles 
were run at one time and restarted 2 times for a total of 60 design cy- 
cles. The slower convergence of the design process with the increase 
in the number of design variables is attributable to the poor perfor- 
mance of quasi-Newton methods for aerodynamic design problems 
with many design variables. 5 As seen in Fig. 8, the desired pressure 
distribution is obtained reasonably well. However, slight waviness is 
noted in the final pressure distribution. The need for curvature con- 
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straints on the geometry is apparent. The final grid is shown in Fig. 9. 



Figure 7. Initial grid for case 2. 



Figure 8. Initial and final pressure distributions for case 2. 



Figure 9. Final grid for case 2. 

Summary and Concluding Remarks 

A two-dimensional design optimization methodology is de- 
scribed. This research is an extension of the work in Ref. 1 to in- 
clude a turbulence model for viscous flows. However, a discrete ad- 
joint approach is used instead of the continuous adjoint approach so 
that the sensitivity derivatives are more consistent with the flow 
solver. The turbulence model is strongly coupled with the flow equa- 
tions, and the accuracy of the derivatives is demonstrated through a 
comparison with derivatives obtained by Finite differences. A few 
examples are presented to demonstrate the methodology. 

In this regard, several recommendations are offered. First, the 
slow convergence of the second test case, in which 75 design vari- 
ables were used, shows that the quasi-Newton method is insufficient 
for problems with many design variables because a large number of 
design iterations is required before a good approximation of the Hes- 
sian can be obtained. Even then, with many design variables, this 
Hessian may remain inaccurate because much of the information is 
obtained much earlier in the design process and may not represent 
the Hessian in the vicinity of the minimum. However, direct compu- 
tation of the Hessian for turbulent Navier-Stokes design cases is not 
currently very efficient or practical because it requires the solution 
of a linear system of equations for each design variable as well as 
one for the adjoint. (See e.g., Refs. 33, and 40.) Methods that ap- 
proximate the Hessian, such as described in Ref. 5, should be thor- 
oughly evaluated and extended to viscous flows. Other methods, 
such as pseudo-time techniques, 36 have been demonstrated for in vis- 
cid flow computations 19 and should be examined for applicability to 
viscous computations as well. In addition, the technique employed 
in Refs. 21, and 22 should also be further evaluated. This technique 
is essentially a time-marching technique in which the gradients are 
smoothed at each step. For two-dimensional flows, this technique is 
similar in application to that of the preconditioning method de- 
scribed in Ref. 5. However, for three dimensions, the technique in 
Ref. 5 requires the solution of an extra field equation. Finally, 
Ta’asan has shown in Ref. 37 that designing for the slopes of the ge- 
ometry instead of the location of the surface presents a design that is 
easier and faster to converge. The use of slopes and curvatures in- 
stead of points as design variables should, therefore, be considered. 
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